#!/bin/bash
#SBATCH --job-name=simulate43
#SBATCH --partition=cpu
#SBATCH -N 1
#SBATCH --ntasks-per-node=8
#SBATCH --mail-type=all
#SBATCH --mail-user=None
#SBATCH --output=/lustre/home/acct-medfzx/medfzx-lkw/project/simulation/log/simulate43%j.stdout.log
#SBATCH --error=/lustre/home/acct-medfzx/medfzx-lkw/project/simulation/log/simulate43%j.stderr.log

cd /lustre/home/acct-medfzx/medfzx-lkw/project/simulation
mkdir -p /lustre/home/acct-medfzx/medfzx-lkw/project/simulation/log
ref=/lustre/home/acct-medfzx/medfzx-lkw/project/CAH/data/refseqMy/GCA_000001405.15_GRCh38_no_alt_analysis_set.fa
outputPath=result/simulate2
BED=/lustre/home/acct-medfzx/medfzx-lkw/project/CAH/data/BED/cah_noname.bed
for fastqDir in data/raw/fastq/22*
do  
    name=$(basename ${fastqDir})
    fastq1=${fastqDir}/${name}.bwa.read1.fastq.gz
    fastq2=${fastqDir}/${name}.bwa.read2.fastq.gz
    # 1 preprocess
    mkdir -p ${outputPath}/fastp/${name}
    fastp -i $fastq1 -o ${fastqDir}/${name}.fastp.read1.fastq.gz \
        -I $fastq2 -O ${fastqDir}/${name}.fastp.read2.fastq.gz  \
        -q 15 \
        -u 40 \
        -e 0 \
        -n 5 \
        -l 50 \
        -p \
        -P 20 \
        -w $(nproc) \
        -g \
        -x \
        -h ${outputPath}/fastp/${name}/${name}.html -j ${outputPath}/fastp/${name}/${name}.json
    module load fastqc/0.12.1
    mkdir -p ${outputPath}/fastqc/${name}
    fastqc -f fastq -o ${outputPath}/fastqc/${name} ${fastqDir}/*.fastq.gz 
    module purge 

    # 2 alignment
    source activate biotools 
    mkdir -p ${outputPath}/bam/${name}
    bamSavePath=${outputPath}/bam/${name}
    bwa mem -t $(nproc) $ref -R "@RG\tID:{name}\tLB:{GRCh38}\tPL:Illumina\tPU:{name}\tSM:{name}" ${fastqDir}/${name}.fastp.read1.fastq.gz ${fastqDir}/${name}.fastp.read2.fastq.gz > ${bamSavePath}/${name}.sam 
    samtools view -@ $(nproc) -bS ${bamSavePath}/${name}.sam  | samtools sort -@ $(nproc) -o ${bamSavePath}/${name}.sorted.bam -
    samtools index -@ $(nproc) ${bamSavePath}/${name}.sorted.bam
    rm ${bamSavePath}/${name}.sam             
    
    # 3 variant calling
    mkdir -p ${outputPath}/vcf/${name}
    vcf=${outputPath}/vcf/${name}/${name}.vcf.gz
    time singularity run ~/singularity/deepvariant.simg \
    /opt/deepvariant/bin/run_deepvariant \
    --model_type WES \
    --ref ${ref} \
    --reads ${bamSavePath}/${name}.sorted.bam \
    --output_vcf ${vcf}  \
    --output_gvcf ${outputPath}/vcf/${name}/${name}.g.vcf.gz \
    --num_shards $(nproc) \
    --regions ${BED} \
    --sample_name ${name} \
    --make_examples_extra_args="min_mapping_quality=1,keep_legacy_allele_counter_behavior=true,normalize_reads=true" \
    # --sample_name

    # # 4 annotation ################ refGene,clinvar_20240611,exac03 #################
    mkdir -p ${outputPath}/vcf/${name}
    annovarOutPut=${outputPath}/vcf/${name}
    annovar=/lustre/home/acct-medfzx/medfzx-lkw/software/annovar
    ${annovar}/table_annovar.pl ${vcf} ${annovar}/humandb -buildver hg38 -out ${annovarOutPut}/${name}_3 -remove -protocol refGene,clinvar_20250721,exac03 -operation g,f,f -nastring . -vcfinput -polish
    
    # MQ output: ${samplePath}/${sampleName}_MQ.csv
    # script/MQ_full.bash -i ${bamSavePath}/${name}.sorted.bam

    bedPath=/lustre/home/acct-medfzx/medfzx-lkw/project/CAH/data/BED/cah.bed
    samtools mpileup --output-MQ -l $bedPath ${bamSavePath}/${name}.sorted.bam > ${bamSavePath}/${name}_MQ.txt
    awk '$1 ~ /^chr/ {print $1 "\t" $2 "\t" $7}' ${bamSavePath}/${name}_MQ.txt > ${bamSavePath}/${name}_MQ_target.txt
    python /lustre/home/acct-medfzx/medfzx-lkw/project/CAH/script/0_totalFlow/WES_script_annovar_LD_MQ/01_uniQGet.py -i ${bamSavePath}/${name}_MQ_target.txt -o ${bamSavePath}/${name}_MQ.csv
    rm ${bamSavePath}/${name}_MQ.txt ${bamSavePath}/${name}_MQ_target.txt 
    # add_LD_MQ
    python script/02_add_LD_MQ.py -i ${annovarOutPut}/*.txt -m ${bamSavePath}/${name}_MQ.csv
    # pseudogene PDP/DP/radio(PDP/DP)  windows-10 PDP/DP/radio(PDP/DP)
    ## step 1 get coverage
    samtools view -b -@ $(nproc) -L ~/project/CAH/data/BED/CYP21A2_CYP21A1P.bed ${bamSavePath}/${name}.sorted.bam > ${bamSavePath}/${name}.target.bam

    samtools depth -b ~/project/CAH/data/BED/CYP21A2.bed -l 50 ${bamSavePath}/${name}.target.bam -@ $(nproc) > ${bamSavePath}/${name}_CYP21A2.depth.tsv
    samtools depth -b ~/project/CAH/data/BED/CYP21A1P.bed -l 50 ${bamSavePath}/${name}.target.bam -@ $(nproc) > ${bamSavePath}/${name}_CYP21A1P.depth.tsv
    ## step 2 add CYP21A2 .. coverage column
    python script/03_add_CYP21A2_CYP21A1P_cov.py -i ${annovarOutPut}/*_3.csv -T ${bamSavePath}/${name}_CYP21A2.depth.tsv -P ${bamSavePath}/${name}_CYP21A1P.depth.tsv
    conda deactivate 
done
